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Groups of eukaryotic cilia and flagella are capable of coordinating their beating over large scales, 
routinely exhibiting collective dynamics in the form of metachronal waves. The origin of this 
behaviour - possibly influenced by both mechanical interactions and direct biological regulation 
- is poorly understood, in large part due to lack of quantitative experimental studies. Here we 
characterise in detail flagellar coordination on the surface of the multicellular alga Volvox earteri, 
an emerging model organism for flagellar dynamics. Our studies reveal for the first time that the 
average metachronal coordination observed is punctuated by periodic phase defects during which 
synchrony is partial and limited to specific groups of cells. A minimal model of hydrodynamically 
coupled oscillators can reproduce semi-quantitatively the characteristics of the average metachronal 
dynamics, and the emergence of defects. We systematically study the model’s behaviour by 
assessing the effect of changing intrinsic rotor characteristics, including oscillator stiffness and the 
nature of their internal driving force, as well as their geometric properties and spatial arrangement. 
Our results suggest that metachronal coordination follows from deformations in the oscillators’ 
limit cycles induced by hydrodynamic stresses, and that defects result from sufficiently steep local 
biases in the oscillators’ intrinsic frequencies. Additionally, we find that random variations in 
the intrinsic rotor frequencies increase the robustness of the average properties of the emergent 
metachronal waves. 

Keywords: eukaryotic flagella, metachronal waves, microhydrodynamics, synchronization, phase 
defects, colloidal oscillators. 


INTRODUCTION 

The eukaryotic flagellum is one of the most highly con¬ 
served structures in biology, providing locomotion and 
fluid transport through the execution of a periodic mo¬ 
tion. From mucociliary clearance in the human respira¬ 
tory tract [1] , to cerebrospinal flows [2] , and the establish¬ 
ment of left-right asymmetry in mammalian embryos [3] , 
the flows generated by this periodic beating are tightly 
coupled to microscale transport in fundamental biologi¬ 
cal processes. Physically separated pairs of flagella can 
synchronize their beating purely through hydrodynamic 
interactions |4], but cilia and flagella usually fulfil their 
tasks through the concerted action of large groups, dis¬ 
playing a remarkably universal tendency to develop syn¬ 
chronous beating patterns known as metachronal waves 
(MWs), large scale modulations of the beating phase [5]. 
Readily observed on the ciliated surface of protists like 
Opalina and Paramecium |S1[7] and even in synthetic bun¬ 
dles of active microtubules [8], MWs are still a mystery: 
their emergence, their properties, and their biological role 
are poorly understood, mostly due to the difficulty of 
studying them in vivo. 

Early qualitative studies of MWs leHii] suggested a 
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mechanical origin for synchronization, and motivated 
theoretical studies of hydrodynamically coupled fila¬ 
ments driven by different internal engines p!2HT6] . While 
these models show a general tendency towards metachro- 
nism, their complexity prevents a simple understanding 
of the underlying mechanism. Minimal models of hydro¬ 
dynamically coupled self-sustained oscillators, instead, 
have the potential to offer insights into the emergence 
of MWs, since in many cases they admit analytical solu¬ 
tions which provide a direct link between model parame¬ 
ters and details of synchrony. Experiments with rotating 
paddles m, light driven microrotors [18], and colloids 
in optical tweezers [IMB, as well as simulations of ro¬ 
tating helices [22|, or spheres driven along fixed [23H27] 
or flexible |28l |29] trajectories, have shown that under 
certain conditions hydrodynamic interactions alone can 
induce phase-locking of two simple oscillators. This lock¬ 
ing is mediated either through wall-modified flows [25] . 
a variable driving force Eaiio], or mechanical elastic¬ 
ity intrinsic in the system [28l [29]. Of these, the last 
two provide the fastest, and possibly equivalently strong 
m, drive towards synchronization. Meanwhile, experi¬ 
ments on the biflagellate alga Chlamydomonas reihnardtii 
|32]-[34], and on isolated pairs of somatic cells from the 
multicellular alga Volvox carteri [4], lend strong support 
to the idea that synchronization among two flagella hap¬ 
pens indeed through the elastic response of flagella to 
shear stress. Even though two flagella can synchronize 
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their motions through direct hydrodynamic interactions 
[4], this pair synchronization does not necessarily guar¬ 
antee group synchronization, nor the emergence of MWs 
[28] . Owing to its large size and ease of visualisation, 
the colonial alga Volvox carteri is an ideal model organ¬ 
ism for the study of flagella-driven flows [35]. Volvox 
comprises thousands of biflagellate somatic cells, embed¬ 
ded within a spherical extracellular matrix shell, beating 
their flagella towards the colony’s posterior with only a 
small deviation 15°) from meridian lines [36]. Here 
we show that these cells display a complex dynamical 
behaviour, with colony-wide MWs propagating towards 
the posterior interrupted by characteristically shaped re¬ 
current phase defects, not previously observed in exper¬ 
iments. For realistic parameters, we show that a mini¬ 
mal model of hydro dynamically coupled identical rotors 
predicts MWs of wavelength similar to those observed 
in experiments. These MWs emerge in a manner es¬ 
sentially independent of boundary conditions, and are 
robust against realistic perturbations in the rotors’ prop¬ 
erties. Including an intrinsic frequency bias derived from 
experiments, the model develops periodic defects simi¬ 
lar to those we report for Volvox. These findings suggest 
that the collective dynamics of Volvox flagella result from 
a competition between a drive towards synchronization, 
based on the oscillators’ elastic compliance, and addi¬ 
tional strains imposed by a large scale bias in flagellar 
properties. 
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FIG. 1. Volvox geometry and flows, (a) A Volvox colony 
held by a micropipette showing flagellated somatic cells (small 
dots) and interior daughter colonies (large circles) growing in 
the posterior half. The dashed line indicates the distance at 
which the kymographs for the components of the fluid veloc¬ 
ity, Ur and U0, are measured, (b) The distribution of radii of 
Volvox colonies used in experiments (total n = 60). (c) Mag¬ 
nitude (colour) and direction (vector) of the time-averaged 
flow field obtained using PIV. (d) Radial component of the 
instantaneous fluid velocity, tZr, at various times throughout 
one flagellar beat cycle. The disturbance propagating towards 
the posterior of the colony is clearly visible. 


METACHRONAL WAVES IN VOLVOX CARTERI 

Volvox carteri f. nagariensis (strain EVE) was grown 
axenically in Standard Volvox Medium [371439] bubbled 
with sterile air. The cultures were hosted in a growth 
chamber (Binder, Germany) set to a cycle of 16 h 
light (100 /rEm“^s“^, Eluora, OSRAM) at 28°C and 
8 h dark at 26°C. Individual Volvox colonies within a 
25 X 25 X 5 mm glass observation chamber filled with 
fresh medium were captured, oriented and held in place 
using two glass micropipettes with 100 /im diameter tips 
housed in pipette holders (World Precision Instruments, 
USA) connected to two manual microinjectors (Sutter 
Instruments Co., USA). Motorised micromanipulators 
(Patchstar, Scientifica, USA) and custom-made rotation 
stages allowed the pipettes to be moved in three di¬ 
mensions and rotated freely along their centrelines. We 
aligned the axis of each colony along the focal plane of 
the 40X Plan Eluor lens (NA 0.6) of a Nikon TE2000- 
U inverted microscope, and recorded 30 s long movies of 
the surrounding flow with a high-speed camera (Eastcam 
SA3, Photron, USA) at 500 fps under bright field illu¬ 
mination. A long pass interference filter with a 10 nm 
transition ramp centred at 620 nm (Knight Optical, UK) 
prevented phototactic responses of the Volvox colonies. 

The projection of the velocity field u onto the focal 
plane was visualised by seeding the fluid with 0.5 jam 
polystyrene microspheres (Invitrogen, USA) at 2 x 10“^ 


volume fraction, and measured using an open source Par¬ 
ticle Image Velocimetry (PIV) tool for Matlab (MatPIV) 
(Fig.0). It was then decomposed into radial and tangen¬ 
tial components n(r, 6>,t) = i4^(r, 6>, t)e^ uo{r,0,t)eo^ 
where the local coordinates are defined with respect to 
the centre of the Volvox colony as shown in Eig. [^a). 
A total of 60 Volvox colonies were selected at random 
from the stock supply, and varied in radius R from 48 to 
251 /rm (mean 144 ± 43/rm), with a distribution shown 
in Fig.Qb). 

The fluid flow time-averaged over the whole dura¬ 
tion to of the movie, (n(r, 6>))^ = V u{r,0R) dt (see 
Eig. [^c)), is described accurately by the modes of the 
‘squirmer’ expansion of the flow field around a sphere @0^ 
Eg modified to take into account the net force exerted 
by the holding pipette. The results are in agreement 
with previous measurements on freely swimming colonies 
[43] . and support the hypothesis that flagellar dynam¬ 
ics is the same in held and freely swimming colonies, a 
fact already established for the closely related unicellular 
species Chlamydomonas [34] . Subtracting the time aver¬ 
age from the instantaneous flow field, u = u — {u)^, high¬ 
lights the flow’s dependence on the local flagellar phase 
within the beating cycle. It is this phase that we wish 
to measure. Eigure shows representative kymographs, 
through the same time interval, for the radial and tan¬ 
gential components of n at r* = 1.3 R. This value of 
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(a) radial velocity Ur xlOO fim/s 




FIG. 2. Properties of metachronal waves. (a) Radial 
Ur{r*^0^t) and (b) tangential ue(r*^O^t) components of the 
flow field, measured at r* = 1.3 R. Phase defects are evi¬ 
denced by white circles, (c) Correlation function C{A0, At) 
for a single representative Volvox colony, (d) its fitted correla- 
tion function cos(fc A6» - cc At) 

and the (e) corresponding error. The scale bar for (c-e) is 
the same, (f) Power spectrum of the autocorrelation function 
C(0, At) (for one Volvox colony), calculated for three distinct 
values of 0. (g) Average beat frequency as a function of polar 
angle 0 for n = 60 different colonies (black) as well as the 
ensemble average (white dotted). 


r*, which corresponds to the dashed circle in Fig. [^a), 
has been found empirically to maximise the kymographs’ 
signals in most experiments. 

Temporal variations in the flow are different for the two 
components (Figs. [^a,b)): Ur exhibits a wave travelling 
from the anterior to the posterior poles of the colony 
(see also Fig. Bd)); uq is closer to a standing wave. 
Such discrepancy should in fact be expected for a MW 
propagating along a spherical distribution of flagella as 
found in Volvox: uq is dominated by the flow induced 
by the flagella in the equatorial region, which are (ap¬ 
proximately) mutually synchronized and are much more 
numerous than those at other latitudes in the colony. 
Peaks in uq would then happen in the middle of the 
power stroke of the equatorial flagella, and should cor¬ 
respond to minima in the magnitude of Ur around the 
equatorial region: this is borne out by the experimental 
kymographs. We confirmed this hypothesis both by vi¬ 
sual inspection of flagellar behaviour and by simulations 
of arrays of model oscillators (see below). The simula¬ 
tions show clearly that only the component of the veloc¬ 
ity perpendicular to the no-slip surface (for us Ur) tracks 
closely the local phase of the oscillators, provided that the 
phase profile changes sufficiently slowly, as is the case for 


Volvox. This finding is also consistent with the fact that 
the flow induced by a force normal to a no-slip plane de¬ 
cays faster - and is then more localised - than that from 
an equivalent tangential force (r“^ vs. r~^ at constant 
distance from the surface). Ur{r*,0^t) can then be used 
as a proxy for the local flagellar phase along the beat¬ 
ing cycle, eliminating the need to track the motion of 
individual flagella on the colony’s surface. Figure ia) 
corresponds to a MW which propagates in the direction 
of the flagellar power stroke, called symplectic, and is 
representative of all the 60 colonies examined (see full ex¬ 
perimental dataset for additional results). Direct inspec¬ 
tion of several colonies at different orientations confirmed 
that Volvox MWs have only a minimal lateral (diaplectic) 
component. The average properties of Volvox MWs can 
be quantified [HI EH] through the normalised autocor¬ 
relation function of Ur{r*,0,t) (Fig. 0c)), C{A0,At) = 
{ur{r'',0 + AO, t -b At) Ur{r'',0, t)) / {^(r*, 6>, t)^) where 
0 denotes averaging over all t and 0. The 
phenomenological functional form CmiAO, At) = 
exp (—^) exp (—^) cos{kA0 — uj At) fits the experi¬ 
mental autocorrelation function remarkably well, with 
average deviation of 0.033±0.011 between the two, across 
all colonies observed. From these fits we estimate the 
mean beating frequency f = uj/27r = 1/T = 32.3±3.1 Hz; 
the autocorrelation decay time r jT = 3.4 ± 2.0, and de¬ 
cay length L = 0.35 ±0.10; and the mean wavenumber 
within the population k = 4.7 ± 0.9 (twice the number 
of complete waves along a meridian), k is positive in 
all our experiments, indicating the ubiquitous presence 
of a symplectic MW which emerges despite the absence, 
in this species, of any direct intercellular connection be¬ 
tween somatic cells. 

By averaging MW properties, we necessarily mask brief 
inhomogeneities in the beating dynamics. These are pri¬ 
marily not random, but rather appear as recurrent de¬ 
fects in the MW pattern, as circled in Fig. 0a). To our 
knowledge, such defects have never previously been ob¬ 
served experimentally, although a similar behaviour was 
reported in recent simulations m- Reminiscent of so- 
called frequency plateaus observed in strips of coupled 
oscillators im, the defects amount to global slip events, 
where the phase difference between two groups of oscil¬ 
lators suddenly increases by a full cycle Hisnaii. In 
the majority of our experiments, the defects are directly 
recognisable from the kymograph (see Fig. ia)). They 
consistently appear around the equatorial region, and 
correspond always to posterior flagella slipping ahead of 
anterior ones. This results in loss of temporal and spatial 
correlation of the MW, and it is responsible for the sur¬ 
prisingly short correlation time (r/T) and length (L) ob¬ 
served. The presence of defects results in the 6>-dependent 
average beating frequency seen in Fig. 0g), which could 
reflect a combination of factors including a physiological 
bias in the intrinsic flagellar beating frequency across the 
colony, or the effect of lower drag due to higher concentra¬ 
tion of somatic cells in the colony’s posterior. Inspection 
of the power spectrum of the autocorrelation function of 
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the kymograph signal at individual values of 6 >, reveals 
that this transition in average beating frequency has a 
characteristic structure. Close to the top or to the bot¬ 
tom of the colony, all cells are phase-locked and beat 
with a single frequency, noticeably different for the two 
groups (the main peaks in the red and blue curves in 
Fig. Hf); the peaks’ width is caused by both intrinsic 
beating variability and measurement error). The transi¬ 
tion region, where defects appear, interpolates between 
these two groups by combining their frequencies. The 
result is a double peak in the power spectrum (e.g. green 
double peak in Fig.j^f)), where the relative height of the 
peaks depends on the proximity to one or the other of the 
phase-locked groups. This phenomenology is character¬ 
istic of all the colonies displaying clearly defined phase 
defects. 

Notice that the frequency bias would by itself gener¬ 
ate a wave propagating in the direction opposite to the 
experimental one, and hence cannot be the fundamental 
mechanism selecting the direction of the observed MW. 
In what follows we discuss the behaviour of a minimal 
model of interacting rotors which spontaneously gener¬ 
ates a symplectic MW, and the effect on these dynamics 
of a large scale bias as found in Fig. [^g). 



FIG. 3. Modelling flagella, (a) Tip trajectory over 10 beats 
of a flagellum of a Volvox somatic cell, (b) Rotor as a model 
flagellum: a sphere of radius a elastically bound to a circular 
trajectory of radius ro (dashed red) perpendicular to a no¬ 
slip plane, driven by a tangential force in the ^-direction, (c) 
Intrinsic trajectory of an isolated rotor. Steady state limit 
cycle of the sphere above the no-slip wall at z = 0 , for A = 
0.1,1,2, 5,10, oo. Perturbations from the circular trajectory 
have been magnified by a factor of 100 so that the shape 
is clearly visible. The evolution of the (d) radius and (e) 
geometric phase through one beating period are also shown. 


HYDRODYNAMIC MODEL FOR INTERACTING 
FLAGELLA 

The surface-mounted somatic cells of Volvox are a few 
tens of microns apart, and their flagella are therefore 
more nearly in the weak-coupling limit than the cilia of 
historically studied organisms such as Paramecium. It 
is thus appropriate to model the fluid disturbance pro¬ 
duced by their operation as a multipole expansion m, 
of which we will only keep the mode with the slowest 
spatial decay - the Stokeslet - representing the effect of 
a point force. This flow is analogous to the far field of 
a rigid sphere pulled through the fluid. It has recently 
been shown that representing a Volvox flagellum as a sin¬ 
gle Stokeslet provides an accurate representation of its 
flow field down to distances of ^ 10 /im, smaller than 
those typically separating cells within colonies (^ 20 fim) 
[4]. High-speed tracking of the flagellar waveform com¬ 
bined with resistive force theory also confirms that the 
distributed forces associated with the motion can be well 
represented by a single point force which periodically tra¬ 
verses a closed loop [4]. Inspired by the trajectories of 
flagellar tips in Volvox^ and following an approach similar 
to others [?7ll29] . a beating flagellum will thus be mod¬ 
elled as a small sphere of radius a elastically bound to a 
circular trajectory of radius ro by radial and transversal 
springs of stiffnesses A and r] respectively, and driven by 
a tangential force of magnitude (see Fig. The 

position of the sphere is given hy x = P s^,r, 0), 
where s = (r sin(0), r cos(0)). The prescribed trajec¬ 
tory, defined by ((" = 0 ,r = ro), is perpendicular to a 
no-slip plane at z = 0 representing the surface of Volvox^ 


and its centre x^ is at a distance d from the plane. The 
proximity of the no-slip boundary causes an asymmetry 
in the sphere’s motion which induces a net flow, thus 
mimicking power and recovery strokes of real flagella. 

The sphere’s velocity v follows the force balance re¬ 
quirement of Stokes flow: 7 (x) • v = —A(r — ro) — 
vC Here 7 = ■y{x) = 70 [I + ( 9 a/ 162 :) (I + 

e,e,)+0((a/zf}] is the friction tensor associated with 
motion near the no-slip wall [46], and 70 = Girfia is the 
drag on a sphere of radius a in an unbounded fluid of 
viscosity /i. Because the drag is not isotropic, the tra¬ 
jectory of an isolated sphere will deviate slightly from 
the prescribed circle {Sr/vQ ^ 0.5% here) to a new limit 
cycle (r(t), 0(t), C(t)) of period T, which we obtain di¬ 
rectly from the simulations for each set of parameters. 
This limit cycle will represent the unperturbed state 
of the single rotor. Figure shows the polar coordi¬ 
nates (r(t), 0 (t)) of such limit cycles, for various val¬ 
ues of the dimensionless ratio A = This 

parameter governs the deformability of the rotors’ or¬ 
bit: larger (smaller) values of A mean stiffer (softer) 
radial recoil compared to the driving force, and result 
in smaller (larger) deviations from r = tq. For Volvox 
flagella, with bending rigidity n = 4 x N m^ [28| 

and length L ^ 10 /im, small deflections of the fila¬ 
ment would be resisted with an effective spring constant 
A n/L^ = 4 X 10“^Nm“^. Approximating the driv¬ 
ing force as = 27rro7o/T, and taking a ^ 1 /im, 

T ^ 1/33 s and d/vQ ^ 1 [29|, the normalised spring con¬ 
stant is estimated as A ^ 0.1. 
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From the limit cycle, we define the phase of the rotor, 
$(0), such that ^ = cj = 27r/T [47]. The presence of the 
wall implies that the phase 4> is different from the angle (j) 
even for completely rigid orbits. The phase will be useful 
when studying the dynamics of coupled rotors, since any 
perturbations from constant evolution of 4> will be solely 
due to hydrodynamic interactions between different ro¬ 
tors. For a system of N rotors, the net hydrodynamic 
force acting on the rotor is 


Fi = -7i 


• ^ ^ (j{xjiXi) ■ Fj 


( 1 ) 


where G(ccj,cc^) is the Green’s function that couples 
the different spheres in the presence of the no-slip wall 
j48l [49] . The time evolution of the system of N moving 
spheres is then given by 


r R 
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Ti 
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y’drive 

-A (rj - ro) 
-vCj 


( 2 ) 

The 3N X 3N matrices F, M and R are defined in terms 
of their constitutive 3x3 blocks Fij = 6ij'j{xi), = 
5ijl^ (l — 6ij)‘^{xi) • G(ccj, cci), and R^j = where 

the columns of the rotation R^ are given by the vectors 


60 ., Br- and bq^ respectively [29] . 


PAIR SYNCHRONIZATION 

Pair synchronization is the building block of large scale 
coordination, and it will be studied here for two identi¬ 
cal rotors spaced a distance I apart, at an angle [3 in 
the x-y plane (Fig. |^a)). For sufficiently stiff springs, 
hydrodynamic interactions induce only small perturba¬ 
tions around the rotors’ limit cycles, and the whole sys¬ 
tem can then be described just in terms of the two 
phases ( 4 >i, 4 > 2 ). Combining the phases into their sum, 
Tj = 4>2 + 4>i, and difference, A = 4>2 — 4>i, provides 
a convenient description of the synchronization state of 
the system, which is simply characterised by A(t). To 
develop intuition, we work almost exclusively within this 
phase-oscillator regime, although simulations use the full 
dynamical system of Eq. (H). 

Figures |^b-e) show the evolution of two coupled 
spheres oriented so that one is directly downstream of 
the other = 7 r/ 2 ), for four different values of A. In 
every case, the rotors are able to achieve a phase-locked 
state by perturbing one another from their limit cycles, 
and the system converges to a steady state in which 
A = A^, independently of the initial condition. The 
fact that A^ < 0 means that the rotor downstream lags 
behind the one upstream, a condition which would be 
naturally expected to produce a symplectic metachronal 
wave along a strip of rotors. We will see below that this 
intuition is not necessarily correct. 
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FIG. 4. Synchronization dynamics of pairs of rotors, (a) 
Schematic diagram of a pair of interacting rotors, each of 
which has three spatial degrees of freedom; geometric phase 
trajectory radius Vi and out-of-plane displacement (not 
shown). The rotor-rotor spacing I and planar angle ^ are de¬ 
picted. (b-e) Simulations of interacting pairs of rotors, for 
various spring stiffness A = Xd/ In each case, the sys¬ 
tem converges to a steady state phase difference with A < 0, 
during a time interval that increases with A. The parameters 
used are a/d = 0.01, ro/d = 0.5, /3 = 7r/2, l/d = 2. 


While As depends very little on the precise value of 
A, this is not the case for the timescale for phase-locking, 
which is seen to increase for progressively stiffer rotors. 
To understand the origin of this dependence, we write the 
evolution of the phase difference A = uj~ -f T)(A, E) and 
sum E = -h 5'(A, E) in terms of effective couplings D 
and S', where = UJ 2 A uji. When the phase difference 
evolves much more slowly than the sum, the equations 
can be averaged over the characteristic timescale of E to 
obtain the time-averaged system: 

A = cj-+5(A) (3) 

E = a;+ -h S(A) (4) 

The crossing of curves in Fig.j^b) indicates that for that 
value of A the evolution of the system is not uniquely de- 
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FIG. 5. Phase oscillator reduction, (a) Coupling functions i^(A) and aS'(A) for a pair of interacting spheres. The corresponding 
parameters are a/d = 0.01, ro/d = 0.5, /3 = 7r/2, l/d — 2^ K — Xd/ = 1. (b) Amplitude of coupling functions for different 
spring stiffnesses. Computed amplitudes of (c) D{A)/uj and (d) S{A)/uj, as well as (e) the stable equilibrium phase difference 
As, shown as functions of rotor-rotor spacing 2 < l/d < 35 and offset angle 0 < yd < 7r/2. 


termined by A: for these parameters the time-averaged 
equations do not provide a good approximation of the dy¬ 
namics. Conversely, each panel in Figs.|^c-e) contains a 
set of curves that, from a given value of A, approach the 
asymptote in essentially the same manner, supporting 
the system’s description in terms of the single variable A. 
Figure l^a) shows the coupling functions D{A) and S{A) 
determined numerically from simulations for one set of 
parameters, and nondimensionalised by the average beat 
frequency u {u = ui = uj 2 in this case). These curves are 
characteristic of the whole parameter range we have stud¬ 
ied ( 0.1 < A < 10). To great accuracy they are, respec¬ 
tively, —Do sin(A —A^) and So cos(A —A^). The system 
has a stable fixed point at A = A^ < 0 , which also cor¬ 
responds to the maximum in S{A): when phase-locked, 
the rotors increase their average beat frequency owing to 
a lower effective drag [ 22 l [28] . Apart from an increase in 
speed, this cooperative drag reduction does not influence 
the synchronization of a single pair, but it will have a 
major impact on the collective state of larger groups of 
rotors. Figure [^b) shows the dependence of the ampli¬ 
tudes Do and So on the effective stiffness A. While So 
is independent of A, reflecting the fact that drag reduc¬ 
tion is a purely hydrodynamic effect, the amplitude 
which governs the convergence towards the stable fixed 
point, scales as Do r\j A ^ (Fig. I^b)). This is reflected 
in the A-dependence of the rate of convergence towards 
phase-locking seen in Fig. [^b). Stiff {Do < So) and soft 
{Do > So) regimes can be recognised, and the position 
of the system with respect to this divide can be tuned 
simply by adjusting the parameter A, notably by chang¬ 
ing the radial stiffness of the orbits. In the asymptotic 
limit as A ^ oo, phase-locking does not occur (for these 
parameters) - a phenomenon represented by Do ^ 0 . As 
reported above, we also find that A 5 is independent of 
A. 

The two most significant parameters governing the spa¬ 
tial configuration of the pair of rotors are the rotor-rotor 


spacing I and the angle P (see Fig. Ba)). Figures [^c-d) 
show the dependence of the coupling strengths {Do, So) 
on {l,P)'. both decrease ^ with increasing /, while 
at fixed separation they increase monotonically as (3 goes 
from 0 (side-by-side) to 7 r /2 (aligned downstream). This 
behaviour follows semi-quantitatively from the (average) 
magnitude of the flow generated by a single rotor, which 
in the far field is proportional to s\n{l3)/l^ [49|. As, 
which is consistently negative, shows instead a much 
weaker dependence across almost all of the region of pa¬ 
rameter space explored (see Fig. |^e)). It decays slowly 
with increasing I {^ for //d > 10 and (3 ^ 7 r/ 2 ) and 
is fairly insensitive to changes in /3, except for a nar¬ 
row region around ^ showing strong dependence on 
the position of the two rotors {As must be an odd func¬ 
tion of p). The phase lag at synchronization is for the 
most part almost unaffected by the detailed configuration 
of the two rotors, in sharp contrast to previous models 
[25l [28] . Simulations for a range of other physical param¬ 
eters, including sphere size and trajectory radius, show 
qualitatively identical results. 

Influence of a variable driving force. — Through¬ 
out its beating cycle, it is plausible that a flagellum will 
exert a changing force on the surrounding fluid and this 
variation, as opposed to waveform compliance, has been 
proposed as a major factor leading to phase locking m- 
Indeed, force modulation can by itself establish synchro¬ 
nization in pairs of model oscillators with rigid trajec¬ 
tories m, and recent studies m have shown that in 
specific circumstances this can dominate over synchro¬ 
nization mediated by orbit compliance. Here we con¬ 
sider the effect of a modulated driving force = 

/o(l + C sin(z/0 + 0o)) in onr system. The average am¬ 
plitude fo is used to define the nondimensional stiffness 
Ao = Xd/fo. Figure]^ shows the coupling functions D{A) 
and S{A) for C = 0.5, u G { 0 , 1 , 2 }, and (fo ^ { 0 , 7 r/ 2 }. 
Notice that v = 2 corresponds to a class of functional 
forms highly effective in establishing synchronization in 
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FIG. 6. Nondimensionalised coupling functions D{/S)/uj (solid) and S{/S)/uj (dotted) for a pair of interacting spheres subject 
to the driving force ~ 1 + O.5sin(z/0 + 0o) for various values (z/, 0o)- Plots are shown with (a) Ao = 1 and (b) Ao = 50. 

Detuning the rotors: evolution of the phase difference A for (c) A = 0.1 and (d) A = 1. Results are shown for various values 
of Q = 2uj~ . Positive and negative values of D are shown in green and red respectively, while the black curves represent 
Q = 0 (identical rotors). All figures correspond to the parameters a/d — 0.01, ro/d = 0.5, /3 = 7r/2, l/d = 2. 


circular rotors [271131]. For Ao = 1 (and smaller), modu¬ 
lating the driving force has no significant bearing on the 
results. While changes in Dq can be significant (< 30%; 
but Do > So always), the difference is mainly far from the 
phase-locked state, and both A^ and So do not vary ap¬ 
preciably from their values at constant driving force. In 
this regime synchronization is dominated by orbit compli¬ 
ance. The corresponding coupling functions for Aq = 50, 
approximating the rigid case, are presented in Fig. [^b). 
So is essentially unchanged, but the phase-locking dy¬ 
namics is severely altered. Both and jA^j have been 
drastically reduced. The most robust synchronization 
appears for (z/, 0o) = (2, 0), with a 75% reduction in Do, 
and this is accompanied by a lag of just jA^j ^ 0.03 (see 
Fig.j^b) inset). As we will see, a system with these char¬ 
acteristics does not generate a MW in a finite chain of 
rotors, an essential requirement for any realistic model 
of metachronism. Of course, there may be specific fixed 
trajectory shapes and force profiles that give rise to an 
appropriate coupling. However, orbit compliance is capa¬ 
ble of facilitating rapid phase-locking in a manner essen¬ 
tially independent of the driving force, in particular for 
Ao ^ 0 . 1 . Additional simulations were undertaken with 
varying rotor geometries, including changes in /3 and /, 
but the results are qualitatively unchanged. 

Detuning the beating frequencies. —Even within 
the same group of cells, different flagella will generally 
have slightly different frequencies. Figure [^g) shows that 
this is the case for Volvox colonies, with a ^ 10% spread 
in (apparent) beating frequencies. Here we study the 
impact of such inhomogeneity on synchronization of a 
pair of rotors. We focus on frequency differences caused 
by variations in driving force. The difference between the 
rotors will be characterised by the fractional frequency 
difference ft = 2uj~ = (co ’2 where O = {uJi + 

Co’ 2)/2 is the mean frequency. Figure |^c,d) shows the 
evolution of the phase difference A for a pair of rotors, 
beginning at A = 0 for t = 0 , for various values of ft. 


The curves have been computed for A = 0.1 and 1 , where 
A is defined in terms of the average driving force. The 
parameters corresponding to the smallest A provide the 
closest approximation to real Volvox cells within this type 
of model. 

Phase-locking is clearly possible for in a finite inter- 
vai around zero, aibeit at different vaiues of A^. Fig¬ 
ure [^c) demonstrates that stabie phase-iocking can oc¬ 
cur for A = 0.1 even when the rotors possess a |11| ^ 6 % 
difference in their intrinsic beat frequencies, which com¬ 
pares weii with the measured vaiue of 10 % for fiageiia 
isoiated from Volvox [4] However, for the iarger vaiue of 
A = 1 , the phase-iocking is iess robust, and the system 
can support a maximum detuning of oniy ^ 0.7% before 
drifting indefiniteiy. As D{A) —Dq sin(A — A^), the 
overaii evoiution of the phase difference can then be writ¬ 
ten as A = ft — Do sin(A — A^) (time is rescaied using 
the average beating frequency). Stationary points occur 
for As (ft) = A^ + arcsin( 0 /Do), which has a soiution 
oniy for \ft\ < Do. Because Do ^ A“^, phase-iocking wiii 
occur beiow a threshoid mismatch |ll|max ^ A“^. This 
estimate, based on the time-averaged phase dynamics, is 
supported by the fuii hydrodynamic simuiations, which 
show that an order of magnitude increase in A (from 0.1 
to 1 ) corresponds to an order of magnitude decrease in 
the threshoid detuning (from ^ 6 % to ^0.7%). 


GROUP SYNCHRONIZATION 

Constant driving force. —Aithough necessary, pair 
synchronization is not a sufficient condition to guarantee 
that a group of rotors wiii generate MWs. This fact is ev¬ 
idenced by studying the behaviour of a finite linear chain 
of rotors, a configuration in which any reaiistic modei 
of metachronai coordination shouid be abie to produce 
robust MWs. We describe here the behaviour of iinear 
chains of 30 identicai rotors, aiigned in the downstream 
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FIG. 7. Metachronal wave development, (a) The phase along the chain of rotors and (b) the neighbouring phase difference 
^i+i — are shown as functions of time. The symplectic MW is established within 15 beats for this value of A = 0.1. (c) 
Phase profile of an array of 30 spheres after t/T = 1200 beats. Results are shown for A = Xd /= 0.1, 0.5,1,1.5,2. All 
results correspond to parameters a/d = 0.01, ro/d = 0.5, Ijd = 2. 


direction, with {ajd^r^^jd^ljd) — (0.01,0.5,2). These 
are meant to mimic the behaviour of cells along a single 
meridian of Volvox^ if they could be observed in isolation. 
The effective stiffness A will be varied from 0.1 {Volvox- 
like) to 2 . Similar studies were performed for chains of 
length ranging from 3 to 50 rotors, and 2 < I/d < 20: the 
qualitative features of the dynamics remain the same. 

Figure shows that for A = 0.1 a steady symplectic 
MW develops from random initial conditions within ^15 
beats. As in a single pair, this phase-locking timescale is 
stiffness dependent, and within our range of A it varies 
from tens to hundreds of beating cycles, independently of 
initial conditions, eventually diverging for A ^ oo. Fig¬ 
ure [^c) shows the phase profiles along the chains after 
1200 beats, representative of the general behaviour. We 
see that the symplectic MW, defined as V i, 

morphs into a chevron pattern as the rotors’ stiffness 
is increased. For this configuration, A 1 marks the 
boundary between the two different patterns. Phase lock¬ 
ing into a chevron has been previously reported for simi¬ 
lar linear systems misni which were not able, however, 
to generate also a stable metachronal coordination. 

Qualitatively, chevron formation arises from the com¬ 
petition between drag reduction at synchronization and 
the tendency of neighbouring rotors to phase-lock at a 
given 7 ^ 0. Drag reduction is more pronounced for 
rotors in the chain’s interior, which will tend to lead those 
at the extremes. This would result in a large fraction of 
the chain having a phase profile opposite to the one re¬ 
quired by metachronal coordination. In the stiff regime 
{Do < So)^ drag reduction dominates, and leads to a 
chevron-like steady state phase profile. However, for suf¬ 
ficiently compliant trajectories {Dq > So) the drive to¬ 
wards locking at a given phase lag is strong enough to 
overcome the effect of drag reduction, and nearest neigh¬ 
bours will synchronize with lags of the same sign as A^. 
In the present system this results in a symplectic MW. 
In fact, even within the MW regime, long range hydro- 
dynamic interactions between rotors at all separations I 
will contribute to determine the steady-state phase lags 
between nearest neighbours. Because As ^ and 


Do ^ we can generally expect the magnitude of 
this phase lag to be smaller than for an isolated pair 
at the equivalent separation. For example, while a pair 
with A = 0.1 would phase-lock with A^ = —0.47 (other 
parameters as in Fig. [^, the average nearest neighbour 
phase lag that is actually realised in a linear array is sig¬ 
nificantly less, A^ = —0.17. The nonlocal interactions 
are clearly important in determining the overall system 
dynamics. With these parameters, an array of 30 ro¬ 
tors should provide the best approximation to a single 
meridian of Volvox cells from a real colony. Our simula¬ 
tions predict that the system should develop a symplectic 
MW, indeed observed experimentally, with a wavenum¬ 
ber k 2. Given the simplicity of the model, we believe 
that this is in good agreement with the experimental av¬ 
erage value oik = 4.7 ± 0.9. 

Inhomogeneous driving forces and defects.— 
The experimental results in Fig. |^g) demonstrate a 
spatially-dependent intrinsic frequency bias among the 
flagella in Volvox. The average shape of this bias is 
captured very well by the functional form ///av = 1 + 
a[tanh( 6 (^ — c)) — tanh( 6 ((i — c))], chosen so that the av¬ 
erage frequency in ^ G [0, tt] is independent of the am¬ 
plitude a. For a chain of 30 rotors, we utilise the fit¬ 
ted driving force /f™//o = 1 + C[tanh(0.1684 x i — 
1.6813) —0.3604] Te^. The parameter fo is chosen to cen¬ 
tre the distribution of resulting effective spring constants 
Ai = around A = 0.1. The amplitude C is a pa¬ 

rameter controlling the extent of the frequency bias. The 
value €i = €i {a) is a random number chosen for each rotor 
from a normal distribution with mean of zero and stan¬ 
dard deviation a, and represents polydispersity in the 
driving forces along the chain. The appropriate parame¬ 
ters for Volvox^ Cy = 0.0724 and ay = 0.0279, were fit¬ 
ted directly from Fig.j^g). The dynamics of rotor chains 
are explored with 0 < C < 0.1 and a G {0,0.0279}, the 
results of which are shown in Fig. 

For sufficiently small frequency biases {C < 0.007), the 
system still converges to a phase-locked state, and the 
only effect is to distort slightly the MW (see Fig. [^a)). 
This shows that the MW is robust, a result further con- 
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FIG. 8. Inhomogeneous driving forces and the appearance of defects. Rotor chains possess a Vo/i’ox-inspired intrinsic driving 
force = 1 + C X (tanh(0.1684 x i — 1.6813) — 0.3604) + a where Ci = u{cr) is a random number chosen from a normal 

distribution with mean of zero and standard deviation a. The frequency bias along the chain is governed by the parameter 
C. (a-d) Kymographs showing the phase, sin(4>i), and phase difference, sin(4>i+i — 4>i), for various values of C, with cr = 0 
(no polydispersity). An example phase defect is circled in white, (e) Kymographs for C — 0.02 with a — 0.0279. (f) Intrinsic 
frequency distributions (black) for 60 different realisations with C = 0.02 with a = 0.0279, along with their average (yellow). 
The magenta curve corresponds to the frequency profile used in (e) and the inset (g) shows power spectra of sin(4>d for 
oscillators 1 (red), 8-9 (green) and 30 (blue) in this chain, (h) Extracted metachronal wavenumber k. Simulations without 
polydispersity {a = 0, black, 19 simulations) and with polydispersity (cr = 0.0279, red, 210 simulations) are shown. Other 
parameters are given hy a/d — 0.01, r^jd — 0.5, Ijd — 2 and centred at A = 0.1. 


firmed by simulations presented in the supplementary 
materials, which explore independently the effects of un¬ 
biased polydispersity and a linear frequency bias. For 
C > 0.007, the effect of the bias is the emergence of de¬ 
fects concentrated around the region of maximal slope in 
the frequency profile (see Figs.|^b-d)). These defects are 
situated in the central/upper part of the chain - the same 
region as for the experimental results in Fig. [^a) - but 
exhibit somewhat slower dynamics than the experimen¬ 
tal defects, lasting approximately 5 to 6 beats, compared 
to 3 to 4 for the experiments. The defects separate two 
regions of symplectic metachronal coordination, a small 
one in the upstream portion of the chain and a larger one 
in the downstream half. In particular, it is evident that 
the bottom third of the strip, corresponding to the region 
of smallest bias, is always coordinated along a symplec¬ 
tic MW. This reflects our experimental observations (see 
Fig. 2 (a)), and it is true for C < 0 . 02 . Above that value, 
this region tends to shrink, going to just 6 oscillators 
for C = 0.08. For values of C large enough to induce 
defects but <0.03 the power spectrum of individual os¬ 
cillators shows the same qualitative behaviour observed 


in the experiments: two well defined frequencies for the 
top and bottom groups of phase-locked oscillators, and 
in the middle a transition region characterised by power 
spectra with a double peak. The relative power carried 
by the two peaks depends on how close the oscillator is to 
one or the other of the phase-locked groups. The transi¬ 
tion region is where the defects appear, as is the case for 
the experiments as well. Altogether, these results sug¬ 
gest that C = 0.02 is an appropriate choice to mimic 
experimental behaviour in the model. 

The effect of including polydispersity in the driving 
forces was also investigated, using the measured value 
from Volvox {a = ay = 0.0279) in the simulations. For 
various values of C, sets of simulations were conducted, 
each a different realisation of this random polydisper¬ 
sity (210 in total, across 9 values of C). Figure [^f) 
shows the intrinsic frequency profile for 60 different sim¬ 
ulations with C = 0.02 (black) together with their aver¬ 
age (yellow), and the underlying bias (white). The ma¬ 
genta curve corresponds to the kymographs presented in 
Fig. I^e). Even in the presence of polydispersity, the 
system displays a symplectic MW on average. 
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FIG. 9. Two-dimensional lattice, (a) Photograph of a patch of flagella on the surface of Volvox (Goldstein). Phase evolution 
for a two-dimensional square grid of 50 rotors with (b) A = Xd /= 2 and (c) the Volvox-\ike parameter A = 0.1. In each 
case, the phase is shown relative to that of the top, central rotor. Other parameters are: orbit radius ro/d = 0.5, rotor-rotor 
spacing Ijd = 2 and sphere size ajd — 0.01. Simulations in b-c have been conducted with free boundary conditions, (d) Average 
lateral and (e) streamwise phase profiles for free boundary conditions (green dashed) and periodic boundary conditions in the 
x-direction (red dashed). 


The correlation function used to characterise the mean 
experimental wave properties is also utilised here to mea¬ 
sure the average metachronal wavenumber in the numer¬ 
ical simulations. For each simulation, the wavenumber k 
was extracted, the results of which are shown in Fig.j^h). 
Without polydispersity (a = 0), increasing the frequency 
bias C results in a reversal of the metachronal wave, with 
k ^ —2 for C > 0.03. However, including random vari¬ 
ations in the intrinsic rotor frequencies (a = 0.0279) in¬ 
creases the wavenumber significantly. For C = 0.03, the 
metachronal wave remains symplectic, with wavenumber 
k 1. Incorporating polydispersity in this fashion in¬ 
creases the frequency of phase defects, which allow relax¬ 
ation of the rotor chain to the locally-preferred symplec¬ 
tic state. The repeated disruption of this metachronal 
wave due to polydispersity in the rotor properties allows 
the wave to sustain clear symplectic coordination on av¬ 
erage, even in presence of a frequency bias that would 
otherwise induce a strong reversal of wave direction. 

Two-dimensional array of rotors. — Figure |^a) 
shows a patch of Volvox cells on the colony’s surface. 
Although the distribution is not entirely regular, the in¬ 
tercellular spacing is almost uniform. Following this ob¬ 
servation, we will investigate in this section the dynamics 
of a two-dimensional array of 10 x 5 identical rotors {x 
and y axes respectively), arranged on a square lattice of 
constant spacing l/d = 2. The rotor orbits are parallel 
to X, with vojd = 0.5, ajd = 0.01, and a radial stiffness 
0.1 < A < 2 . 

Commencing from random initial phases, the system 
with A = 2 moves over several hundred beats towards 
a phase profile in which the interior rotors lead (see 


Fig. [^b)). This is analogous to the profile exhibited for 
the linear chain in Fig.j^c), whereby drag reduction of in¬ 
terior rotors leads to a chevron phase profile. Figure [^c) 
shows that for A = 0.1 the system develops a symplec¬ 
tic MW within a few tens of cycles with average phase 
difference of = —0.19 between adjacent rotors with 
the same ^-coordinate. This is very similar to the value 
found for a linear chain {As = —0.17), indicating that 
the one-dimensional configuration is able to capture the 
salient features of the system. Figures |^d-e) show the 
average lateral and streamwise phase profiles respectively 
(green dashed curves). 

The same behaviour was observed when imposing pe¬ 
riodic boundary conditions along y, mimicking the peri¬ 
odic arrangement of surface-mounted flagella on Volvox. 
Starting from random initial conditions, all our simula¬ 
tions for A = 0.1 converged rapidly to a symplectic MW 
with rotors along lines of constant x locking in-phase. 
There was almost no measurable effect on the streamwise 
metachronal wave profile (red dashed curve in Fig. ie)), 
though the slight phase variations in the ^-direction were 
suppressed. This periodic arrangement can also sustain 
a discrete set of lateral MW components. Our simula¬ 
tions show that a MW with a total phase accumulation 
of ±27r along the periodic direction {y) is locally stable, 
and higher windings ±27rn should also be possible, al¬ 
though their domain of stability is likely to shrink with 
increasing integer n. These states did not develop in 
any of our simulations starting with random initial con¬ 
ditions, so the probability to realise them over the purely 
symplectic wave is likely to be very small. 
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CONCLUSIONS 


This paper presents a detailed experimental analysis of 
the global flagellar dynamics on the surface of the multi¬ 
cellular green alga Volvox carteri. As a model organism, 
Volvox has the fundamental advantage that its flagella 
are well separated from each other and hence much more 
in the low-coupling limit than for any other organism 
where metachronal coordination has been previously re¬ 
ported (e.g. protozoa, ciliated epithelia). This property 
allowed us to compare successfully flagellar dynamics in 
Volvox with a minimal model of hydro dynamically cou¬ 
pled elastic rotors. The model develops spontaneously 
a symplectic MW of a wavenumber that compares well 
with the one observed experimentally. This MW is robust 
against inhomogeneities in the rotors’ intrinsic frequen¬ 
cies, variations in their spatial arrangement, and choice 
of boundary conditions. These are all essential properties 
for a plausible model of metachronal coordination. We 
have shown the existence of MWs which are consistently 
symplectic on average, but display at the same time char¬ 
acteristic recurrent defects. These might be connected 


to the observed bias in the distribution of beating fre¬ 
quencies across the colony, which by itself would favour 
a wave propagating in the opposite direction. We tested 
this hypothesis within our minimal model, and observed 
the emergence of recurrent phase defects in the same lo¬ 
cation and with the same power spectral signature as 
those observed experimentally. Altogether, these results 
show that different properties of the underlying oscilla¬ 
tors can combine and give rise to much more complex 
metachronal wave dynamics than previously assumed. 
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Supplementary Material 

The results of various simulations are shown here 
for chains of rotors above a no-slip wall, with inhomo¬ 
geneities in their driving forces. We explore the effects 
of introducing a frequency bias of the same functional 
form as the measured frequency distribution in Volvox 
(see Fig. 2(g)), and also additional polydispersity in the 
rotor driving forces. This frequency profile in Fig. 2(g) 
was fitted and for a chain of 30 rotors, the appropriate 
driving force is given by = 1 + C[tanh(0.1684 x 

i - 1.6813) - 0.3604] with Cy = 0.0724 and 1 < i < 30 
is the rotor index. The parameter /o is chosen to cen¬ 
tre the distribution of resulting effective spring constants 
Ai = around A = 0.1. Figure 2 (g) also shows 

that the frequency profiles for individual Volvox colonies 
are distributed around the mean curve. The standard de¬ 
viation of this spread was measured to be ay = 0.0279 for 
Volvox. In the following sections, various driving forces 
and degrees of polydispersity are considered. Multiple re¬ 
alisations of each parameter set are shown, highlighting 
the general behaviour of each configuration. 


INCREASING FREQUENCY BIAS 

(0.001 <C < 0.1) 

NO POLYDISPERSITY {a = 0) 

For small values of C, the shape of the MW is per¬ 
turbed, yet remains phase-locked and symplectic in na¬ 
ture. For C > 0.007 phase defects begin to emerge 
among the interior rotors. These defects are periodic 
for 0.007 < C < 0 . 02 , but for larger values of C, defects 
emerge at various points in the chain and interact to pro¬ 
duce more complex phase dynamics. In every simulation, 
there exists a region of oscillators at each end which are 
locally phase-locked. 



FIG. 10. Chains of rotors with an intrinsic frequency bias. 
The driving force along each chain is given by = 

1 -h C[tanh(0.1684 x i — 1.6813) — 0.3604] for various values 
of C G [0.001,0.1]. The value for Volvox is given by Cy = 
0.0724. Other parameters are given by A = 0.1, a/d = 0.01, 
ro/d = 0.5, l/d — 2. The phase profile and phase difference 
from each simulation are shown, as well as the intrinsic rotor 
frequency profile (red curves) and average measured frequency 
profile from the simulations (blue curves). As the frequency 
bias is increased, defects in the middle of the chain emerge, 
becoming more abundant for larger values of C. 









































































14 


FIXED FREQUENCY BIAS (C = 0.02) 
POLYDISPERSITY (a = 0.0279) 


FIXED FREQUENCY BIAS (C = 0.0724) 
POLYDISPERSITY (a = 0.0279) 



FIG. 11. Various realisations of polydispersity. The driv¬ 
ing force along each chain is given by = 1 + 

0.02[tanh(0.1684 x i — 1.6813) — 0.3604] + ei where Ci = ei(cr) 
is a random number chosen from a normal distribution with 
mean of zero and standard deviation a = av — 0.0279. Other 
parameters are given by A = 0.1, a/d = 0.01, r^jd — 0.5, 
lld — 2. The phase profile and phase difference in each simu¬ 
lation are shown. The colour scale includes —1 (blue) through 
to +1 (red). The results corresponding to 10 of the 60 simu¬ 
lations are shown here, and are representative of the general 
behaviour observed. 


FIG. 12. Various realisations of polydispersity. The driv¬ 
ing force along each chain is given by = 1 + 

C[tanh(0.1684 x i — 1.6813) — 0.3604] + ei where C = Cv = 
0.0724 and Ci = ez(cr) is a random number chosen from a nor¬ 
mal distribution with mean of zero and standard deviation 
a = av = 0.0279. Other parameters are given by A = 0.1, 
ajd — 0.01, rojd = 0.5, Ijd = 2. The phase profile and phase 
difference in each simulation are shown. The colour scale 
includes —1 (blue) through to +1 (red). The results corre¬ 
sponding to 10 of the 60 simulations are shown here, and are 
representative of the general behaviour observed. 





































































































































































NO FREQUENCY BIAS (C = 0) 
POLYDISPERSITY (a = 0.0279) 


Correlation parameters 
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FIG. 13. Various realisations of polydispersity. The driv¬ 
ing force along each chain is given by = 1 + 

where Ci = is a random number chosen from a nor¬ 

mal distribution with mean of zero and standard deviation 
cr = cry = 0.0279. Other parameters are given by A = 0.1, 
ajd — 0.01, rojd — 0.5, Ijd = 2. The phase profile and phase 
difference in each simulation are shown. The colour scale 
includes —1 (blue) through to +1 (red). The results corre¬ 
sponding to 10 of the 30 simulations are shown here, and are 
representative of the general behaviour observed. Symplectic 
MWs persist on average, even in the presence of polydisper¬ 
sity. 


In this section, the correlation parameters are ex¬ 
tracted for chains of rotors with varying degrees of poly¬ 
dispersity. The driving force for the rotor is taken to 
be /f"™//o = (1 + Cj) where €i = €i{a), as in Fig. 


but for various values of a. The parameter /o is cho¬ 
sen to centre the distribution of resulting effective spring 
constants = Xd /around A = 0.1. With these 
parameters, the threshold detuning for an isolated pair 
is |Q|max — 6 %. Here we want to investigate the extent 
to which this applies to MWs, and extract the correlation 
parameters for the emergent metachronal waves. 

Altogether, we performed 150 simulations of arrays of 
30 rotors, each run for approximately 3000 cycles start¬ 
ing with random initial conditions. For each simulation, 
the limit cycle of individual rotors was independently de¬ 
termined, and used to define their phases and beating 
periods The normalised standard deviation 

of the set of periods = std({T^}?£;|^)/Tav characterises 
the variation in the set of beating periods for each simu¬ 
lation. Here Tav is the average of {T^}. As 'ip is increased, 
the MW that would develop in a system of identical ro¬ 
tors is perturbed initially only slightly {ip = 0.009; rotors 
globally phase-locked), then more heavily {ip = 0.028; 
rotors locally phase-locked), until eventually the MW is 
realised only on average and any locking is only local 
{ip = 0.052). Kymographs were used to estimate the av¬ 
erage MW properties (T^r^L^k) for each simulation, as 
described previously for the experiments. 

Figures [T^ a-b) show that, generally, both the autocor¬ 
relation length and time, L and r, decrease significantly 
as ip is increased. In these plots, the systems with r 
larger than the duration of the simulation are coloured 
green, while the remaining points are shown in blue. The 
autocorrelation time plot reveals that these two groups 
cluster around markedly distinct values of r. Simulations 
in the first group have r/T ^ 10^ — 10^ and seemingly 
independent to the value of ip. They develop a steady 
MW despite the finite dispersity, but only up to ip ^ 3%. 
This value corresponds to the threshold detuning we en¬ 
countered for two rotors, because for a pair ip = Cl{2. In 
the second group, r decreases with ip. Figure p^c) dis¬ 
plays the spread of wavenumbers k. The mean wavenum¬ 
ber {k = 2 . 1 , dashed line) compares well with the value 
for identical rotors. With a finite dispersity in the driv¬ 
ing forces, the wavenumbers spread in a manner similar 
to the distribution observed experimentally for Volvox 
(Fig. p^d)). A few systems at large values of ip dis¬ 
played negative wavenumbers (^ 6 % of the simulations), 
but overall the simulations show that the present model 
robustly exhibits symplectic MWs, even with variations 
in rotor properties. 
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FIG. 14. (a) Correlation decay length L and (b) time r/T 
scaled by period are shown as functions of the spread in 
beating periods b- (c) Histogram of extracted wavenumbers. 
Values of b are the same as the preceding two plots. Re¬ 
sults correspond to 150 simulations of V = 30 spheres with 
a/d = 0.01, ro/d = 0.5, l/d = 2 and centred at A = 0.1. 
(d) Histogram of wavenumbers in Volvox colonies [29] . 


LINEAR FREQUENCY BIAS 
NO POLYDISPERSITY {a = 0) 


To this point, the Volvox-inspired intrinsic frequency 
profile has been used, with the parameter C controlling 
the extent of this bias along the chain. Chains of rotors 
with a linear bias in their intrinsic frequencies are stud¬ 
ied here. The rotors are actuated with a driving force of 
the form /f™//o = 1 + T) [i — (A + l)/ 2 ] for various 
values of D G [—0.0003,+0.0005], so that their intrin¬ 
sic beat frequency varies monotonically along the chain. 
Figure [Tb] shows the steady state phase profiles exhibited 
for a number of values of D. This frequency bias can ei¬ 
ther enhance {D < 0) or reduce {D > 0) the slope of the 
metachronal wave, while still permitting convergence to 


a steady state. The direction of the metachronal wave is 
robust, even when the constituent rotors possess a mild 
opposing frequency bias 1 % difference between the 
end rotors). 

There are significant qualitative differences between 
this system, and the one presented in Fig.[^ For the lin¬ 
ear profile here, a very small overall frequency difference 
between the end rotors {D = +0.0005) results in sup¬ 
pression of the symplectic metachronal wave. Effective 
hydrodynamic coupling for each rotor does not extend to 
more than a few neighbours away, and so the mild linear 
frequency profile does not provide “weak points” in the 
coupling. Conversely, the intrinsic frequency profiles in 
Fig. vary nonlinearly. Defects emerge in the middle 
of the rotor chain where the neighbouring frequency dif¬ 
ference is largest, allowing relaxation to the symplectic 
metachronal wave. 



FIG. 15. Steady state metachronal wave for chains of rotors 
with an intrinsic frequency bias. Inset shows the intrinsic 
frequency of rotors along each chain. The functional form 
of the driving force is linear, and given by fo = 1 + 

D[i — {N for various values of D G [—0.0003,+0.0005]. 
Other parameters are given by A = 0.1, a/d = 0.01, ro/d = 
0.5, l/d = 2. 




























